For generating most, but not all figures in the manuscript.
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
[30m── [1mAttaching packages[22m ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──[39m
[30m[32m✓[30m [34mggplot2[30m 3.3.0 [32m✓[30m [34mpurrr [30m 0.3.3
[32m✓[30m [34mtibble [30m 3.0.0 [32m✓[30m [34mdplyr [30m 0.8.5
[32m✓[30m [34mtidyr [30m 1.0.2 [32m✓[30m [34mstringr[30m 1.4.0
[32m✓[30m [34mreadr [30m 1.3.1 [32m✓[30m [34mforcats[30m 0.5.0[39m
[30m── [1mConflicts[22m ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:cowplot’:
stamp
The following objects are masked from ‘package:dplyr’:
intersect, setdiff, union
The following objects are masked from ‘package:base’:
date, intersect, setdiff, union
library(mgcv)
Loading required package: nlme
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
source("UVP_2017_library.R")
Loading required package: future
Parsed with column specification:
cols(
Cruise = [31mcol_character()[39m,
Station = [31mcol_character()[39m,
`mon/dd/yyyy` = [31mcol_character()[39m,
`hh:mm` = [34mcol_time(format = "")[39m,
`Longitude [degrees east]` = [32mcol_double()[39m,
`Latitude [degrees north]` = [32mcol_double()[39m,
`Bottom Depth [m]` = [32mcol_double()[39m,
`Pressure [db]` = [32mcol_double()[39m,
`Temperature [degrees C]` = [32mcol_double()[39m,
`Temperature 2 [degrees C]` = [32mcol_double()[39m,
`Salinity [psu]` = [32mcol_double()[39m,
`Salinity 2 [psu]` = [32mcol_double()[39m,
`Fluorescense [mg/m^3]` = [32mcol_double()[39m,
`Beam Transmission [%]` = [32mcol_double()[39m,
PAR = [32mcol_double()[39m,
`Oxygen [umol/kg]` = [32mcol_double()[39m,
`Oxygen [% saturation]` = [32mcol_double()[39m
)
theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
These read back in in UTC, since write_csv saves everthing into that time zone.
options(readr.default_locale=readr::locale(tz="Mexico/General"))
bes<- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
bds <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
ues <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m,
depth = [32mcol_double()[39m,
psd_gam = [32mcol_double()[39m,
vol = [32mcol_double()[39m,
sizeclass = [31mcol_character()[39m,
lb = [32mcol_double()[39m,
ub = [32mcol_double()[39m,
binsize = [32mcol_double()[39m,
TotalParticles = [32mcol_double()[39m,
nparticles = [32mcol_double()[39m,
n_nparticles = [32mcol_double()[39m,
biovolume = [32mcol_double()[39m,
speed = [32mcol_double()[39m,
flux = [32mcol_double()[39m,
flux_fit = [32mcol_double()[39m,
GamPredictTP = [32mcol_double()[39m
)
uds <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
Specify the base of the photic zone (which is inside of the OMZ) and the bae of the OMZ
PhoticBase <- 160
OMZBase <- 850
library(scales)
Attaching package: ‘scales’
The following object is masked from ‘package:purrr’:
discard
The following object is masked from ‘package:readr’:
col_factor
#https://stackoverflow.com/questions/30179442/plotting-minor-breaks-on-a-log-scale-with-ggplot
log_breaks = function(maj, radix=10) {
function(x) {
minx = floor(min(logb(x,radix), na.rm=T)) - 1
maxx = ceiling(max(logb(x,radix), na.rm=T)) + 1
n_major = maxx - minx + 1
major_breaks = seq(minx, maxx, by=1)
if (maj) {
breaks = major_breaks
} else {
steps = logb(1:(radix-1),radix)
breaks = rep(steps, times=n_major) +
rep(major_breaks, each=radix-1)
}
radix^breaks
}
}
scale_x_log_eng = function(..., radix=10) {
scale_x_continuous(...,
trans=log_trans(radix),
breaks=log_breaks(TRUE, radix),
minor_breaks=log_breaks(FALSE, radix))
}
#theme_set(theme_bw)
PlotPSDmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) +
scale_shape_manual(name = "Day", values = c(21:25)) +
scale_fill_gradientn(name = "Hour",breaks = c(0, 6, 12, 18, 24), labels = c("0", "6", "12", "18", "24"), colors = c("black", "blue", "white", "orange", "black")) +
labs(y = "Depth (m)", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue")
#theme_set(theme_cowplot)
PlotParticlesmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
scale_x_log10(breaks = c(10, 100, 1000), minor = c(5, 50, 500)) +
#theme(legend.position = "none") +
#scale_x_log_eng()+
labs(y = "Depth (m)", x = "Particles / L") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue") +
theme(legend.position = "none")
plot_grid(
PlotParticlesmany,
PlotPSDmany,
rel_widths = c(1, 2)
)
Removed 266 rows containing missing values (geom_point).Removed 266 rows containing missing values (geom_point).
ggsave("figures/ParticlesPSDMany.png")
Saving 12 x 7.41 in image
ggsave("figures/ParticlesPSDMany.svg")
Saving 12 x 7.41 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
#FSG4 <- gam(tot_nparticles~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.96538 0.09655 92.86 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.0000 1.000 5.663 0.0202 *
s(Day) 1.4473 1.694 2.033 0.0921 .
s(Hour) 0.6402 2.000 0.502 0.2154
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.135 Deviance explained = 17.4%
GCV = 0.68369 Scale est. = 0.64319 n = 69
summary(FSG2)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.96538 0.09727 92.17 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.000 1.000 5.741 0.0194 *
s(Day) 1.469 1.718 2.214 0.0792 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.122 Deviance explained = 15.3%
GCV = 0.68744 Scale est. = 0.65288 n = 69
summary(FSG3)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.9654 0.1004 89.34 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1 1 5.736 0.0194 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.0651 Deviance explained = 7.89%
GCV = 0.71559 Scale est. = 0.69485 n = 69
#summary(FSG4)
summary(FSG1)$r.sq - summary(FSG2)$r.sq # extra R^2 explained by hour
[1] 0.01304448
summary(FSG2)$r.sq - summary(FSG3)$r.sq # extra explained by day
[1] 0.05646094
summary(FSG3)$r.sq # R^2 explained by depth
[1] 0.06511156
There is no statisticlly significant affect of time on particle number. If I take all of the time variables out and just compare to depth, there is a relationship to depth p = 0.02. But the R^2 is only 6.5%. Pretty weak
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG4 <- gam(psd~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.90483 0.01809 -215.8 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.848 1.977 88.789 <2e-16 ***
s(Day) 1.506 1.756 0.953 0.4853
s(Hour) 1.361 2.000 2.260 0.0424 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.716 Deviance explained = 73.6%
GCV = 0.024631 Scale est. = 0.022591 n = 69
#summary(FSG2)
#summary(FSG3)
summary(FSG4)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.90483 0.01814 -215.2 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.845 1.976 88.630 <2e-16 ***
s(Hour) 1.524 2.000 2.559 0.0409 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.714 Deviance explained = 72.8%
GCV = 0.024244 Scale est. = 0.022709 n = 69
summary(FSG3)$r.sq # R^2 of depth
[1] 0.6929817
summary(FSG4)$r.sq - summary(FSG3)$r.sq # Improvement from adding hour of the day
[1] 0.02133449
summary(FSG1)$r.sq - summary(FSG4)$r.sq # Improvement from then adding day of the week
[1] 0.001481684
PSD varies with depth, but doesn’t statistically relate to hor orday. Comparing the R2 values from models tells us that you explain 69% of varience with depth.
PlotNParticlesEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
#geom_path(aes(x = tot_nparticles)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
labs(x = "Particles/L", y = "Depth (m)") +
theme(legend.position = "none") +
scale_shape_manual(values = c(1:5)) +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue")
PlotNParticlesEP
I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)
PlotPSDEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
geom_path(aes(x = psd_gam)) +
geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_color_manual(values = c("gray20", "brown")) +
scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue")
PlotPSDEP
plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))
Removed 1211 rows containing missing values (geom_point).Removed 1211 rows containing missing values (geom_point).Removed 1211 row(s) containing missing values (geom_path).
ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")
Saving 10 x 4 in image
ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")
Saving 10 x 4 in image
Large and spall particle number, flux and size
mainParticleComponents <- bds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
select(project, profile, depth,
tot_nparticles, small_nparticles, big_nparticles,
tot_psd = psd, small_psd, big_psd,
tot_flux_fit, small_flux_fit, big_flux_fit) %>%
pivot_longer(cols = -c("project", "profile", "depth")) %>%
separate(name, c("size", "meas")) %>%
mutate(meas = recode(meas, nparticles = "particles/L")) %>%
mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))
Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].
PlotFlx <- mainParticleComponents %>%
filter(meas != "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
PlotPSD <- mainParticleComponents %>%
filter(meas == "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))
Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).
ggsave("figures/BigVsSmall.svg")
Saving 7.29 x 4.5 in image
ggsave("figures/BigVsSmall.png")
Saving 7.29 x 4.5 in image
Flux small and flux tot track so closely because particle fractal dimension alpha, plus flux fractal dimension, gamma > |psd|. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top.
Example particle size distributions
eg_dataline <- bds %>%
filter(profile == "stn_043", depth == 162.5)
eg_slope = eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)
eg_datablock <- bes %>%
filter(profile == "stn_043", depth == 162.5)
eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)
eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)
EgNNP <- eg_datablock %>%
ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")
EgNP <- eg_datablock %>%
ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")
EgTP <- eg_datablock %>%
ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")
plot_grid(EgNNP, EgTP, labels = c("A", "B"))
Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis
ggsave("figures/ExamplePSD163m.png")
Saving 7.29 x 4.5 in image
ggsave("figures/ExamplePSD163m.svg")
Saving 7.29 x 4.5 in image
Flux attenuation with respect ot depth and time. All extrapolated from the UVP and traps combined.
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_b <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_c <- function(x) {
xout <- gsub("1e", "10^{", format(x),fixed=TRUE)
xout <- gsub("{-0", "{-", xout,fixed=TRUE)
xout <- gsub("{+", "{", xout,fixed=TRUE)
xout <- gsub("{0", "{", xout,fixed=TRUE)
xout <- paste(xout,"}",sep="")
return(parse(text=xout))
}
scale_x_log10nice <- function(name=NULL,omag=seq(-10,20),...) {
breaks10 <- 10^omag
scale_x_log10(breaks=breaks10,labels=scientific_10_c(breaks10),...)
}
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10nice()+
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
#labs(x = "moo", y = "Depth (m)") +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm")) +
geom_segment(aes(y = 160, yend = 160, x = 20, xend = 500), color = "darkgreen", stroke = 0.5)+
geom_segment(aes(y = OMZBase, yend = OMZBase, x = 20, xend = 500), color = "darkblue", stroke = 0.5)#+ geom_hline(yintercept = OMZBase, color = "darkblue")
Ignoring unknown parameters: strokeIgnoring unknown parameters: stroke
pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)
Removed 14 rows containing missing values (geom_point).
pltFlx
#plotly::ggplotly(plt1)
Zooming in on where the action is happening
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse()+
#scale_x_log10() +
scale_x_log10(breaks = c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20)), limits = c(20, 180)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")+
geom_hline(yintercept = 160, color = "darkgreen")
pltFlxZoom
#plotly::ggplotly(plt1)
Rate of change of flux, taken to the fifth root so one can see patterns.
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")+
geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3
#plotly::ggplotly(plt1pos)
Within panel drawing
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","B"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgTop
pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)
I don’t know whats going on below here
pgBottom <- pltDelta3 + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = "white", inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)
Removed 33 rows containing missing values (geom_point).
pgBoth
ggsave("figures/FluxDeepDive.png")
Saving 5 x 9 in image
ggsave("figures/FluxDeepDive.svg")
Saving 5 x 9 in image
Test for day to day and hourly variability in rate of change of flux (fifth root transformed)
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_DayOnly <- gam(pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_NoFifth <- gam(pracma::nthroot(DF/DZ, 1)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05927 -3.102 0.0037 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.6736 1.893 3.885 0.0635 .
s(Day) 1.9071 1.991 5.159 0.0132 *
s(Hour) 0.9147 2.000 0.922 0.1473
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.272 Deviance explained = 35.2%
GCV = 0.16976 Scale est. = 0.14755 n = 42
summary(DFG_DayOnly)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06436 -2.857 0.00682 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Day) 1.879 1.985 4.024 0.0299 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.142 Deviance explained = 18.1%
GCV = 0.18677 Scale est. = 0.17396 n = 42
summary(DFG_NoFifth)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 1) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.028354 0.007893 -3.592 0.000941 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.809e+00 1.963 9.249 0.00138 **
s(Day) 1.855e+00 1.979 3.028 0.05229 .
s(Hour) 1.625e-08 2.000 0.000 0.77687
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.346 Deviance explained = 40.5%
GCV = 0.0029435 Scale est. = 0.0026167 n = 42
# summary(DFG2)
# summary(DFG3)
#
# summary(DFG1)$r.sq - summary(DFG2)$r.sq
# summary(DFG2)$r.sq - summary(DFG3)$r.sq
# summary(DFG3)$r.sq
There is variability with respect to depth, and day and hour Depth p = 0.03 R^2 = 0.088. Add affect of day p = 0.004, extra R^2 = 0.11, Add affect of hour p = 0.02 extra R2 = 0.12
bdsAddTime <- bds %>%
mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)
DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_DayOnly <- gam(pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_NoFifth <- gam(pracma::nthroot(DF/DZ, 1)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05998 -3.065 0.00407 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.672 1.893 3.955 0.0606 .
s(Day) 1.907 1.991 4.661 0.0190 *
s(Hour) 0.774 2.000 0.658 0.1961
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.254 Deviance explained = 33.4%
GCV = 0.1732 Scale est. = 0.15112 n = 42
summary(DFG_DayOnly)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06485 -2.835 0.00722 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Day) 1.873 1.984 3.391 0.0404 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.129 Deviance explained = 16.8%
GCV = 0.18962 Scale est. = 0.17665 n = 42
summary(DFG_NoFifth)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 1) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.028354 0.007924 -3.578 0.000979 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.809e+00 1.963 9.247 0.00138 **
s(Day) 1.827e+00 1.970 3.361 0.06222 .
s(Hour) 1.396e-08 2.000 0.000 0.71927
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.341 Deviance explained = 39.9%
GCV = 0.0029647 Scale est. = 0.0026375 n = 42
# summary(DFG2)
# summary(DFG3)
#
# summary(DFG1)$r.sq - summary(DFG2)$r.sq
# summary(DFG2)$r.sq - summary(DFG3)$r.sq
# summary(DFG3)$r.sq
As above, but this time continuous day and hour
Plot of the gams above
#plot.new()
FluxGamPlot <- function(){
par(mfrow = c(2,2))
plot(DFG1)
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}
FluxGamPlot()
png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)
FluxGamPlot()
dev.off()
png
2
(u mol C / m^3 / day)
disagFig <- bds %>% filter(project == "ETNP") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
#labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day)), y = "Depth (m)") +
labs(x = paste("Deviation from Model", expression((μmol/m^3/day)))) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = PhoticBase, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue")
disagFig
#ggsave("..figures/FluxSizeShift.svg"
ggsave("figures/FluxSizeShift.png")
Saving 6 x 4 in image
ggsave("figures/FluxSizeShift.svg")
Saving 6 x 4 in image
# bdsAddTime <- bds %>%
# mutate(Hour = hour(time), Day = day(time))
bdsAddTime <- bds %>%
mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)
OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(OZG1)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01049 7.665 1.33e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.584 1.827 19.185 4.77e-06 ***
s(Day) 1.771 1.946 8.926 0.00173 **
s(Hour) 1.266 2.000 1.964 0.05140 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.416 Deviance explained = 45.5%
GCV = 0.0082683 Scale est. = 0.0075948 n = 69
summary(OZG2)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01079 7.451 2.85e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.558 1.805 18.393 8.14e-06 ***
s(Day) 1.838 1.974 6.924 0.00415 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.382 Deviance explained = 41.3%
GCV = 0.0085853 Scale est. = 0.0080382 n = 69
summary(OZG3)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01169 6.877 2.59e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.453 1.7 17.31 2.23e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.274 Deviance explained = 29%
GCV = 0.0097831 Scale est. = 0.0094354 n = 69
summary(OZG1)$r.sq - summary(OZG2)$r.sq # Extra from Hour
[1] 0.03410746
summary(OZG2)$r.sq - summary(OZG3)$r.sq # Extrafrom Day
[1] 0.1074734
summary(OZG3)$r.sq # Depth
[1] 0.2741944
Plot of those gams Figure S10
OSMSGamPlot <- function(){
par(mfrow = c(2,2))
plot(OZG1)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(2,1))
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
}
OSMSGamPlot()
png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 8, units = "in", res = 200)
OSMSGamPlot()
dev.off()
png
2
trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")
Parsed with column specification:
cols(
Class = [31mcol_character()[39m,
Depth = [32mcol_double()[39m,
TrapID = [31mcol_character()[39m,
TrapType = [31mcol_character()[39m,
SampleType = [31mcol_character()[39m,
C_flux = [32mcol_double()[39m,
C_flux_umol = [32mcol_double()[39m
)
UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
Flux = [32mcol_double()[39m
)
UVPFluxOE <- read_csv("dataOut/ObservedVsExpectedFlux.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
tn_flux = [32mcol_double()[39m,
profile = [31mcol_character()[39m,
project = [31mcol_character()[39m,
time = [31mcol_character()[39m,
tot_flux2 = [32mcol_double()[39m
)
fluxMS_distilled_toPlot <- trapFlux3 %>%
mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))
Remove traps where mass spec didn’t work correctly 2-17 150 1-12 73m 1-12 148 2-14 100 |(TrapID == “2-17” & Depth == 150)
fluxMS_distilled_toPlot2 <- fluxMS_distilled_toPlot %>%
filter(!((TrapID == "1-12") | (TrapID == "2-14" & Depth == 100)|(TrapID == "2-17" & Depth == 150)))
#fluxMS_distilled_toPlot2
Traps where mass spec didn’t work.
UVPFluxPlot00 <- UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, shape = TrapType, ID = TrapID),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot2) +
geom_line(aes(x = Flux), size = 1, color = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP Estimate")) + # dummy point for the legend
geom_point(aes(x = tot_flux2), size = 3, shape = 21, color = "white", fill = "black", data = UVPFluxOE) +
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") +
#xlab(expression(Flux µmolC/m^2/day)) +
xlab(expression(paste("x axis ", ring(A)^2))) +
xlab(expression(paste("Flux (µ mol C/", m^2, "/day)"))) +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
)
Ignoring unknown aesthetics: ID
# UVPFluxPlot <- UVPFluxPlot00 +
# geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
UVPFluxPlot00
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
TPPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = TotalParticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "TotalParticles Observed (#)")
nnpPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = n_nparticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Binsize and Volume Normalized Particles (#/L/mm)")
FitPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = nnp_smooth, xmin = nnp_lower, xmax = nnp_upper, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Smoothed - Normalized Particles (#/L/mm)") + geom_errorbar(width = 10, alpha = 0.5)
npLegend <- get_legend(FitPlot + theme(legend.box.margin = margin(0, 0, 40, 200)) + labs(col = expression(log[e](Size (mm)))))
Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
plot_grid(
TPPlot + theme(legend.position = "none"),
nnpPlot + theme(legend.position = "none"),
npLegend ,
FitPlot + theme(legend.position = "none")
)
Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
ggsave("figures/AllParticleSizes.svg")
Saving 10 x 6.18 in image
ggsave("figures/AllParticleSizes.png")
Saving 10 x 6.18 in image
A bunch of upfront smoothing
SameGam <- gam(TotalParticles ~s(log(lb), log(depth), by = factor(time)), offset = log(vol * binsize), family = nb(),
data = bes %>% filter(project == "ETNP", depth <= 2000))
besE <- bes %>% filter(project == "ETNP")
lb_new <- exp(seq(from = log(0.1), to = log(2.1), by = 0.05))
ub_new <- lead(lb_new)
binsize_new <- ub_new - lb_new
lbbs <- tibble(lb = lb_new, ub = ub_new, binsize = binsize_new)
Expanded <- expand_grid(lb = exp(seq(from = log(0.1), to = log(2), by = 0.05)), depth = seq(from = 20, to = 2000, by = 20), time = as.factor(unique(besE$time))) %>% left_join(lbbs, by = "lb")
Pred <- exp(predict(SameGam, Expanded))
ToPlot <- bind_cols(Expanded, nnparticles = Pred) %>% mutate(time = as.character(time)) %>% mutate(nparticles = nnparticles * binsize)
meanBese <- ToPlot %>% filter(lb <= 2) %>% group_by(lb, depth) %>% summarize(nparticles = mean(nparticles), nnparticles = mean(nnparticles))
WBColorMap <- meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = expression(log[10](Particles/m^3/mm))) + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs(y = "Depth (m)", x = "Size (mm)")
WBColorMap
mbGam <- meanBese %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
pWBPSD <- mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue")
pWBPSD
PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- PubDf %>% filter(depth <= 165, depth >= 155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 160, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs( x = "Small particle mass (norm.)")
pWBS
LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 165, depth >=155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 160, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue")
pWBL
WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))
Removed 6 row(s) containing missing values (geom_path).Removed 7 row(s) containing missing values (geom_path).
WBFig5
WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined
ggsave("figures/WBModelValidation.png")
Saving 8 x 6 in image
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxP16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10(limits = c(35, 150),breaks = seq(from = 20, to = 150, by = 20)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
geom_hline(yintercept = 200, color = "darkgreen") +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))
#
#
#
# pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
# pltFlxLegend <- get_legend(pltFlx)
#
pltFlxP16
# #plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3P16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-1, .1), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 200, color = "darkgreen")+
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3P16
#plotly::ggplotly(plt1pos)
osms_p16 <- bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), group = factor(time))) + geom_point(size = 3) + geom_path() + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = "Observed - Modeled Small Particle Flux \n µmol/m^3/day") +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
plotly::ggplotly(osms_p16)
#ggsave("..figures/FluxSizeShift.svg"
plot_grid(
pltFlxP16,
pltDelta3P16,
osms_p16
)
Removed 28 rows containing missing values (geom_point).Removed 28 row(s) containing missing values (geom_path).Removed 32 rows containing missing values (geom_point).Removed 32 row(s) containing missing values (geom_path).Removed 29 rows containing missing values (geom_point).Removed 29 row(s) containing missing values (geom_path).
ggsave("figures/P16FluxRelate.svg")
Saving 8 x 8 in image
ggsave("figures/P16FluxRelate.png")
Saving 8 x 8 in image
Take one profile, attenuate it, and show what it looks like
source("ModelStuff.R")
Parsed with column specification:
cols(
Cruise = [31mcol_character()[39m,
Station = [31mcol_character()[39m,
`mon/dd/yyyy` = [31mcol_character()[39m,
`hh:mm` = [34mcol_time(format = "")[39m,
`Longitude [degrees east]` = [32mcol_double()[39m,
`Latitude [degrees north]` = [32mcol_double()[39m,
`Bottom Depth [m]` = [32mcol_double()[39m,
`Pressure [db]` = [32mcol_double()[39m,
`Temperature [degrees C]` = [32mcol_double()[39m,
`Temperature 2 [degrees C]` = [32mcol_double()[39m,
`Salinity [psu]` = [32mcol_double()[39m,
`Salinity 2 [psu]` = [32mcol_double()[39m,
`Fluorescense [mg/m^3]` = [32mcol_double()[39m,
`Beam Transmission [%]` = [32mcol_double()[39m,
PAR = [32mcol_double()[39m,
`Oxygen [umol/kg]` = [32mcol_double()[39m,
`Oxygen [% saturation]` = [32mcol_double()[39m
)
scan_for_example <- bds %>% filter(project == "ETNP", depth < 500, depth > 200) %>% select(profile, depth, DFP, use_this_DFP, ospsDZ)
#loc_station = "stn_036"
loc_station = "stn_043"
loc_depth = 225
loc_prev_depth = 112.5
allDFPs <- bds %>% filter(profile == loc_station, depth >= loc_prev_depth, depth <= loc_depth) %>% summarize(DFP = prod(DFP), use_this_DFP = prod(use_this_DFP))
loc_DFP <- allDFPs %>% pull(DFP)
loc_use_DFP <- allDFPs %>% pull(use_this_DFP)
for_single_disag <- bes %>% filter(profile == loc_station, depth %in% c(loc_prev_depth, loc_depth)) %>% select(depth, lb, nnp_smooth) %>%
mutate(depth = recode(depth, `112.5` = "Shallow", `225` = "Deep")) %>% # I have no idea how to not hard code this bit
pivot_wider(names_from = depth, values_from = nnp_smooth)
with_disag <- for_single_disag %>%
mutate(Predicted_Deep = remin_smooth_shuffle(Shallow, loc_use_DFP))
#remin_smooth_shuffle(for_single_disag$Shallow,loc_use_DFP)
for_plot_disag <- with_disag %>% pivot_longer(cols = -lb) %>% #filter(lb <= 5) %>%
mutate(name = factor(name, levels = c("Shallow", "Deep", "Predicted_Deep"))) %>%
mutate(name = recode_factor(name, Shallow = "Shallow (112.5m)", Deep = "Deep (225m)", Predicted_Deep = "Predicted Deep (225m)"))
for_plot_disag %>% ggplot(aes(x = lb, y = value, shape = name)) + geom_point() + scale_x_log10() + scale_y_log10() + scale_shape_manual(values = c(1, 6, 3)) + theme(legend.title = element_blank()) + labs(x = "Particle Size (mm)", y = "Normalized Particle Abundance (#/L/mm)")
ggsave("figures/DisagExample.png")
Saving 7.29 x 4.5 in image
ggsave("figures/DisagExample.svg")
Saving 7.29 x 4.5 in image